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Abstract 

The  research  under  this  contract  have  mainly  been  along  two  directions:  (1)  three-dimensional 
(3D)  finite  difference  (FD)  modeling  of  scattering  from  free  surface  topography,  and  (2)  near  real 
time  event  location  using  only  seismogram  envelopes  from  local  networks.  In  (1)  we  have  simulated 
scattering  from  teleseismic  P-waves  using  a  plan  vertical  incident  P-wave  and  real  topography 
from  a  40  x  40  km  area  centered  at  the  NORESS  array  in  south-eastern  Norway.  Snapshots 
and  synthetic  seismograms  of  the  wavefield  show  clear  conversion  from  P  to  Rg  (short  period 
fundamental  mode  Rayleigh)  waves  in  an  area  of  rough  topography  approximately  10  km  east  of 
NORESS.  This  result  is  consistent  with  numerous  observations. 

In  (2)  a  new  method  is  described  for  automatic  epicenter  locations  in  near  real  time  using  short 
period  z-component  data  from  local  seismograph  networks.  The  original  waveform  data  is  bandpass 
filtered,  ’STA  envelope’  parameterized  and  then  resampled  at  a  rate  of  2  Hz.  The  physical  principle 
invoked  for  epicenter  determinations  is  that  of  constructing  the  space-time  image  of  a  source  in 
the  gridded  network  area  on  the  basis  of  P-  and  S- wavelet  intensities.  Since  such  intensities  for  a 
stratified  half  space  is  almost  model  independent  we  did  not  bother  to  use  local  travel  time  tables 
nor  crustal  models.  The  latter  information  is  naturally  needed  for  precise  origin  time  estimates. 
The  method  is  robust  because  and  to  our  surprise,  P-  and  S-intensities  (low  frequency  wavelets) 
vary  smoothly  (ID)  with  distances  at  least  out  to  700  km.  Compared  to  conventional  event  location 
approaches,  we  bypassed  tasks  like  signed  detection,  phase  identification  and  phase  association  - 
in  essence  we  aim  for  joint  event  detection  and  location.  The  method  has  been  very  convincingly 
tested  on  Italian,  German  and  Norwegian  network  data  off-line  but  has  not  been  adapted  to  a  real 
time  operational  environment  nor  attempted  extended  to  teleseismic  recordings.  In  the  former 
case  we  see  no  principal  problems  and  besides  computational  requirements  can  be  handled  by  a 
PC  machine.  The  essence  of  our  automatic  epicenter  location  scheme  is  the  potential  for  fast 
(3-4  min),  effective  and  low-cost  exploitation  of  CTBT  relevant  monitoring  information  from  the 
very  many  local  networks  deployed  in  most  continental  areas.  In  the  future  we  foresee  that  local 
network  Hubs  may  inform  the  IDC  of  a  CTBT  treaty  about  incoming  P-waves  at  teleseismic 
located  stations  and  arrays  from  specific  sources  at  given  locations.  Also,  the  ever  increasing  use 
of  Internet  for  seismic  data  exchange  may  at  least  partly  move  CTBT  seismic  monitoring  into  the 
seismological  community  sphere  and  perhaps  partly  into  the  public  domain. 
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Figure  1:  Leftmost  section  shows  the  topography  in  an  100  x  100  km  area  centered  at  the  NORESS 
array.  The  innermost  40  x  40  km  was  used  in  the  modeling  experiment.  The  dotted  lines  shows 
the  positions  of  the  receivers  in  the  two  profiles  and  the  circle  outlines  the  NORESS  array.  Labels 
are  in  kilometers  and  elevations  are  in  meters  above  mean  sea  level.  The  black  area  is  Lake  Mj0sa 
(123  m  above  sea  level). 

1  Objective 

The  main  goal  of  the  basic  seismological  research  performed  under  this  contract  is  to  improve 
our  understanding  of  wave  propagation  and  scattering  phenomena  for  local  and  regional  distance 
travel  paths.  The  tool  used  here  is  numerical  simulations  through  the  very  general  and  powerful, 
but  computationally  demanding,  finite  difference  technique. 

Another  research  goal  is  aimed  at  exploiting  the  information  potential  of  local  network  record¬ 
ings  in  the  context  of  CTBT  monitoring,  the  first  step  here  is  the  development  of  a  fast,  cost- 
efficient  method  for  automatic  event  locations  using  such  network  data.  To  ensure  robustness, 
conventional  analysis  tasks  like  signal  detection,  phase  identification  and  association  were  avoided. 

2  Basic  research  results  (  B.o.Ruud  and  s.Hesthoim ) 

2.1  3D  FD  modeling  of  scattering  from  topographic  relief 

The  scattering  of  teleseismic  P  wave  energy  into  Rg  have  been  extensively  studied  and  documented 
using  data  from  the  NORESS  array  in  south-eastern  Norway  (Bannister  et  al,  1990;  Gupta  et  al, 
1993;  Hedlin  et  al,  1991,  1994).  For  our  3D  FD  modeling  experiments  we  have  obtained  digital 
elevation  data  for  an  area  of  100  x  100  km  centered  on  the  NORESS  array  (Fig.  1).  Due  to  the 
huge  computer  memory  requirements  of  3D  FD  methods,  we  have  so  far  been  restricted  to  a  model 
of  size  40  x  40  x  35  km  with  0.2  km  sampling. 
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Figure  2:  Snapshots  of  the  vertical  component  of  the  particle  velocity.  The  upper  snapshots 
are  horizontal  sections  at  the  free  surface  and  the  lower  snapshots  are  vertical  sections  along 
a  west-east  profile  through  the  center  of  the  model.  The  time  in  the  lower  left  corner  of  each 
horizontal  section  is  the  time  from  when  the  plane  vertical  incident  P-wave  was  reflected  from  the 
surface.  The  straight  wavefronts  parallel  to  the  model  boundaries  are  artificial  reflections  from 
the  absorbing  boundaries.  These  artificial  reflections  are  seen  also  in  the  vertical  sections.  Note 
the  strong  scattering  source  located  at  Bronkeberget  (2  km  N,  10  km  E)  which  is  also  highly 
prominent  in  real  record  analysis  (Bannister  et  al,  1990). 


The  method  used  to  implement  free  surface  topography  in  the  3D  FD  code  is  an  extension 
of  the  2D  method  described  by  Hestholm  and  Ruud  (1994).  In  all  the  examples  shown  here  the 
incoming  wave  is  a  vertical  incident  plane  P-wave  simulating  a  teleseismic  short  period  P-phase. 
The  center  frequency  of  the  Ricker  wavelet  is  2.5  Hz,  the  P-wave  velocity  of  the  homogeneous 
medium  is  6.0  km/s  and  the  Poisson’s  ratio  is  0.25.  The  absorbing  boundary  condition  used  along 
the  sides  and  bottom  of  the  model  is  the  so  called  ‘exponential  damping’  technique. 

The  main  problem  encountered  in  our  test  runs  was  instabilities  which  started  near  the  surface 
in  areas  of  steep  and  rough  topography.  The  instabilities  appear  as  surface  waves  of  exponentially 
increasing  energy  which  slowly  propagate  out  from  the  starting  points.  The  first  modeling  attempt 
with  unsmoothed  real  topography  gave  very  unstable  results  (infinite  amplitude)  within  a  second 
after  the  wave  front  reached  the  surface.  In  our  next  attempt  the  elevations  were  first  multiplied 
by  0.5  in  order  to  reduce  the  topography.  The  results  were  stable  and  this  is  also  consistent  with 
previous  2D  FD  modeling  results.  Also  in  the  latter  case  an  west-east  profile  running  through 
NORESS  was  unstable  with  real  topography  but  stable  with  half  elevations.  Snapshots  of  the 
wavefield  (vertical  component  of  the  particle  velocity)  are  shown  in  Fig.  2.  A  dominant  feature 
of  the  snapshots  is  the  low-frequency /long- wavelength  artificial  reflections  from  the  absorbing 
boundaries.  Although  the  exponential  damping  technique  is  the  most  efficient  absorbing  boundary 
method  we  have  tested,  the  case  of  an  wave  incident  at  an  angle  of  90°  with  the  boundary  is  the 
most  difficult  one.  Fortunately,  the  artificial  reflections  are  frequency  dependent,  with  longer 
wavelengths  than  most  of  the  surface  waves  scattered  by  the  topography,  and  can  therefore  be 
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Figure  3:  Seismograms  extracted  at  the  free  surface  of  the  model.  Left  section  shows  the  vertical 
component  of  the  particle  velocity  along  a  30  km  long  west-east  profile  through  the  center  of  the 
model.  The  right  section  shows  3-component  seismograms  from  a  receiver  3.6  km  east  of  the 
center. 

removed  by  filtering  and  image  processing  techniques.  As  seen  from  Fig.  2  the  scattered  surface 
waves  appear  to  radiate  out  from  secondary  point  sources  which  coincide  with  areas  of  high 
topographic  gradients  (Fig.  1).  The  dominant  scattering  points  are  along  the  steep  valley  side 
east  of ‘Bronkeberget’  about  10  km  east  of  NORESS  (Bannister  at  al.,  1990).  Also  in  seismograms 
from  the  west-east  profile  the  P-to-Rg  scattering  from  this  area  is  clearly  seen  (Fig.  3). 

2.2  Recommendations  and  future  plans  -  3D  FD  synthetics 

Even  for  small  scale  experiments  3D  FD  synthetic  wavefield  analysis  provide,  as  demonstrated 
here,  improved  insight  and  a  better  understanding  of  surface  scattering  phenomenon.  We  find  it 
particularly  gratifying  that  the  strong  P-to-Rg  scattering  from  the  Bronkeberget  hill,  observed  by 
Bannister  et  al.  (1990)  through  analysis  of  NORESS  recordings,  can  be  realistically  synthesized. 
Another  interesting  feature  is  that  for  certain  spatial  low-frequency  variants  of  the  local  NORESS 
topography  the  Rg  wave  propagation  seems  to  change  abruptly  over  relatively  small  distances. 
Probably,  this  is  due  to  strong  directionality  dependence  of  the  scattering  from  some  topographic 
features.  Such  wavefield  characteristics  are  sometimes  observed  in  NORESS  and  GERESS  record¬ 
ings  -  at  some  sensors  the  Rg-phase  is  prominent  while  hardly  visible  at  nearby  sensors  less  than  a 
kilometer  away.  In  other  words,  also  small  scale  crustal  features  may  contribute  to  blocking  effects 
as  often  observed  for  Rg-  and  Lg-phases.  Many  of  such  wave  field  phenomena  are  out  of  range  for 
2D  FD  synthetic  experiments  hence  our  emphasize  on  continuous  research  efforts  on  3D  synthetic 
simulation  of  crustal  wave  field  propagation. 

Recently,  we  have  modified  the  code  to  run  on  parallel  computers  and  we  are  now  in  the 
testing  stage.  Future  directions  of  our  research  will  be  to  increase  the  model  size  to  also  include 
the  scattering  areas  SE  of  Lake  Mj0sa.  Furthermore,  in  order  to  compare  directly  with  NORESS 
recordings  it  is  necessary  to  allow  for  non- vertical  incident  waves  and  to  use  source  time  functions 
derived  from  teleseismic  P  beams.  Although  the  physical  dimensions  of  3D  models  still  will  be 
quite  small  compared  with  2D  FD  models,  the  use  of  parallel  computer  technology  opens  new 
avenues  to  study  of  3D  scattering  phenomena  in  fields  like  random  media  and  corrugated  interface 
scattering. 
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3  Automatic  event  location  using  local  network  data 

(  G.A.  Ryzhikov,  M.S.  Biryulina  and  E.S.  Husebye  ) 

As  a  preamble  to  deciding  on  a  local  event  location  strategy  we  attempted  to  formulate  sort 
of  ’boundary  conditions’  for  approach.  Firstly,  to  tie  the  analysis  to  the  parameterized  waveforms 
make  sense  also  scientifically  because  signal  waveforms  are  dissimilar  for  station  separations  ex¬ 
ceeding  a  few  kilometers.  Likewise,  an  approach  in  principle  similar  to  conventional  ones  including 
signal  detection,  phase  identification  and  phase  association  was  rated  unattractive  because  corre¬ 
sponding  parameters  extraction  in  an  automatic  mode  is  unlikely  to  be  robust  (Ruud  et  al,  1993). 
Finally,  for  local  events  there  is  no  need  to  ignore  a  priori  information  embedded  in  the  event 
records  namely  the  presence  of  the  two  P-  and  S-wave  modes  and  with  the  latter  being  relatively 
energetic.  Incorporating  amplitude  information  would  naturally  reduce  the  severity  of  the  phase 
association  problem  due  to  fast  signal  amplitude  decay  with  epicenter  distance. 

From  the  above  considerations,  we  formulated  the  following  epicenter  location  strategy: 

1.  Analysis  restricted  to  parameterized  waveform  data  say  tied  to  rms  or  to  the  STA  detector 
parameter  definition  (sampling  rate  0. 5-1.0  sec) 

2.  Extraction  of  ID  (distance  dependence  only)  signal  attributes  to  ensure  robustness  of  anal¬ 
ysis. 

3.  A  grid  search  approach  giving  an  option  for  including  a  priori  information  in  on-line  analysis. 

4.  Epicenter  location  criteria  based  in  minimizing  differences  between  observed  and  expected 
seismic  signal  attribute  features. 

5.  ’Expected’  attributes  to  be  anchored  in  seismic  wave  propagation  theory. 

3.1  Theory  -  event  localization  strategy  realization 

Out  of  the  various  task  listed  above,  that  of  transforming  the  original  waveforms  or  parameterized 
waveform  records  in  such  a  manner  that  distinct  signal  attributes  could  be  extracted  reliably, 
proved  difficult  to  realize  in  practise.  In  particular  the  increasing  scattering  contribution  to  the 
observed  wavefield  with  time  tend  to  obscure  distinct  signal  features  for  larger  (A  >  300  km) 
epicenter  distances.  With  the  failure  of  our  empirical  experiments  in  minimizing  distinct  signal 
attributes,  we  examined  many  works  on  wave  propagation  in  stratified  media  for  a  clue  to  this 
problem  (e.g.,  see  Brekhovskikh  1960,  Aki  and  Richards  1981,  Kennett  1983,  Flatte  et  al,  1979). 
On  this,  presumably  adequate  theoretical  basis,  we  proceeded  in  a  manner  reminiscent  of  Chap. 3 
and  4  in  the  Aki  and  Richards  book.  In  other  words,  the  starting  point  being  the  far  field 
displacement  for  a  time-space  distributed  source  in  an  inhomogeneous  media  expressible  as: 

=  Gtj(x,  x';  t  -  r)si(x',t)  (1) 

where  ipi  is  the  ith  displacement  component  at  location  x  (x  =  {xi,X2,x$}  d=  {r,  z}),  r  being 
the  free  surface  location  and  z  is  depth)  and  time  t;  G,j  is  the  Green  operator  (tensor  function) 
and  Si  is  the  body  force  equivalent  at  position  x'  and  time  r.  The  Green  function  is  complicated, 
depending  essentially  on  an  infinite  number  of  medium  parameters  and  as  such  not  attractive  in 
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a  real  time  analysis  environment.  Hence  we  rewrite  eq.  (1)  for  an  isotropic  stratified  medium  in 
the  following  form: 

<p°  =  G°s  (2) 

with  index  zero  signifying  the  simplified  medium.  We  may  rewrite  the  general  Green  function  in 
terms  of  that  for  the  simplified  media,  namely 

G  =  G°  +  6G  (3) 

where  G°  —  G°(|x  —  x'|;  t  —  r).  Introducing  the  Lame  operator  L  on  G: 

LG  =  6(x  -  x')6(t  -  r)  (4) 

where  the  L-operator  in  time-frequency  representation  is  of  the  form: 

L  =  -pu26ik  -  dj{djkidi]  (5) 

with  p—  density,  6  =  Kronecker  delta,  ui  =  angular  frequency  and  c  =  elasticity  tensor.  Using  the 
L  -operator,  we  have  for  6G  (eq.  3): 

6G  = -G°(L  -  L°)G  =  G°VG  (6) 

where  the  potential  V  =  L°  —  L.  Obviously,  seismic  waveforms  vary  rapidly  between  stations  so 
preferences  are  for  wavefield  intensity  or  amplitude  square  measures,  that  is 

<p<p*  =  GG'ss*  (7) 

where  G*  stands  for  the  operator  being  adjoint  to  G.  Inserting  eqs.(2),  (3),  and  (6)  in  eq.  (7)  we 
get 

=  y»V*  +  G°G°VV*W*  (8) 

where  the  overline  imply  averaging  over  two  separate  recording  positions  and  operator  VV*  is 
supposed  to  be  essentially  nonlocal.  Note  that  equation  (8)  describes  wavefield  intensities  taking 
into  account  that  we  have  omitted  the  expression  6(xi  —  X2)S(ti  —  1 2)  for  equation  (8);  all  terms 
should  be  multiplied  by  the  expression. 

Now,  recalling,  that  our  ’target  zone’  is  the  free  surface  we  are  interested  in  describing  of 
energy  flux  in  lateral  direction,  with  G°  describing  the  propagation  in  a  stratified  halfspace  with 
smooth  perturbations.  V  can  be  modeled  as  an  ensemble  of  random  functions  with  lateral  radius  of 
correlation  being  dominant,  which  gives  the  wave  propagation  in  terms  of  adiabatically  perturbed 
normal  modes.  This  means  that  propagation  in  lateral  direction  r  —  r*  (r*  is  epicenter)  can  be 
treated  in  terms  of  a  parabolic  wave  equation,  and  eq.  (8)  can  be  treated  as  an  equation  of 
evolution  of  intensity  w  =  ipip*6(x  1  —  X2)S(ti  —  t2).  Physically,  we  have  at  this  stage  a  simplified 
expression  for  wavefield  intensity  with  the  term  being  the  energy  flux  in  an  isotropic  medium 

and  the  VV*-  term  giving  diffusion  of  energy  during  propagation.  Note,  the  propagation  of  an 
impulse  signal  in  isotropic  stratified  half  space  has  two  main  wavefield  intensity  components  in 
a  vicinity  of  the  free  surface:  P-primary,  or  ’3D-mode’  and  S-secondary,  or  principal  ’2D-mode’, 
having  specific  group  velocities,  which  are  defined  for  the  P-mode  by  turning  points  at  depth  z  >  0 
and  for  5-mode  is  caused  by  a  turning  point  at  depth  z  =  0  (Kennett,  1983).  Now  the  difference 
Aw  =  w  —  w°  between  intensities  in  an  inhomogeneous  half  space  w  and  in  an  isotropic  stratified 
one  has  to  satisfy  the  diffusion  equation: 

dA w(r,t)  _  2d2Aw(r,t) 

dr  a  dt 2  W 
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where  a  is  the  diffusion  coefficient.  Equation  (9)  gives  the  distance  dependence  of  wavefield 
intensity  for  localized  (point)  source  function. 

Observationally,  seismogram  intensity  may  be  of  form  of  the  STA  detector  (Ruud  et  al,1993), 
and  we  want  to  minimize  the  differences  between  observed  u  and  modeled  w  wavefield  intensities, 
namely  to  find  w  such  that 

M  =  {u—  w\u  —  w)  — ►  min  (10) 

where  u  —  {u( rn,z  =  0;  t)  =  u(rn,t)  =  un\n  =  1,2,..., AT),  where  r„  is  the  n-th  station,  N  is 
number  of  stations  and  w  =  tu(|r  —  r'|,  z  =  0;  t  —  r),  and 

1  N  r 

<“M=f  /  «(»«, 

n= ljAt 


We  may  write  a  model  for  observations  u,-  on  the  base  of  modeled  w  (eq.(8))  in  an  approximate 
manner  by  assuming  linearity  in  terms  of  intensities,  that  is: 

Wn  —  +Swn  =  Wf  +  6wn  (11) 

where  W  =  G°G°*(|rn  —  r'|,  z  =  0;  t  —  r)),  /  =  /(r',t)  is  source  function  squared;  8wn  =  du>(r„,t) 
is  an  error  term,  which  can  be  interpreted  as  a  random  realization  of  Atn  with  covariance  operator 
Cw  =  8wn6wn  =  exp  {a2rndf}  (eq.(9)).  The  equation  =  Wf  is  similar  to  eq.(l)  except  that 
displacement  is  replaced  by  wavefield  intensity  and  at  that  in  a  simplified  manner.  From  eq.(ll) 
we  can  estimate  the  source  intensity  /(r',t): 

/  =  (W*C-1W  +  air1W’C-1u  (12) 


where  a  is  positive  parameter  of  regularization.  Using  a  RT-approximation  (Ryzhikov  and  Troyan 
1992  a,b)  and  presuming  linearity  of  w  with  respect  to  /,  we  may  rewrite  eq.(ll)  for  dimensionless 
values  /  =  ///max  and  u  =  u/umax  in  the  following  form: 


/Or'.  T) 


(  Wfr.r  1  Tu  ) 
{Wv,T  |HV,T)  +  a 


(13) 


where  observed  data  u  are  filtered  by  operator  T  =  exp{— (l/2)a2r„5t2}  (  Fig.  4c)  and  the  unit 
intensity  source  response  Wv\T  =  W(|r  —  r"|,  t  —  t')8(r"  —  r ')8{t'  —  r)  (  an  eq.(l)  analogy).  Albeit 
/  may  be  used  as  an  event  magnitude  measure,  we  prefer  to  work  with  normalized  observational 
data.  This  ensures  that  corresponding  source  intensity  function  /  varies  between  0  to  1.0  and 
plays  the  role  of  correlation  between  the  real  and  virtual  source  intensity  functions.  It  allows  us 
to  introduce  the  source  image  function: 


V’Or',  r)  =  exp  {cr~2[/(r',  r)  -  1]} 


(14) 


where  /( r',  r)  is  the  normalized  value  of  /  in  eq.(13).  The  entropy  of  source  image  contrast  (EnIC), 
defined  by  Biryulina  and  Ryzhikov  (1995)  is  : 

^(r)  =  “  /  Mr'>T))ln/*(r,,r))dr/  (15) 

Jn 

via  a  pseudo-probability  density  function  p(r ',t)  =  (9r'V')2(r',  r)/  Jn(c’r'V’)2dr/.  The  attractive 
feature  of  the  EnIC  is  that  of  beig  a  very  good  measure  for  detecting  a  single  impulse-like  seismic 
event.  A  more  detailed  theoretical  basis  of  this  epicenter  location  scheme  is  given  by  Ryzhikov, 
Biryulina  and  Husebye  (1995,  in  preparation) 
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Figure  4:  Local  event  record  for  the  Italian  station  FVI  at  an  epicenter  distance  of  240  km  (date 
11/02/95  -  Table  1).  a)  Original  waveforms  after  2  -  4  Hz  bandpass  filtering;  sampling  rate  is  50 
Hz.  b)  Parameterized  version  of  the  (a)  waveform  subjected  to  a  simple  non-overlapping  short¬ 
term-average  ’envelope’  transform;  sampling  rate  is  2  Hz.  c)  P-  and  S-intensity  wavelets  derived 
from  the  parameterized  waveforms  shown  in  (b).  The  peak  of  these  wavelets  as  a  function  of 
distance  (Fig.  5)  sure  the  only  information  used  in  our  automatic  location  scheme.  Basic  seismic 
network  traces  used  in  analysis  is  the  (b)-type  parameterized  waveforms  which  are  easily  produced 
in  situ.  Note,  we  have  for  convenience  used  the  classical  P-  primary  and  S-secondary  notations  for 
the  intensity  wavelets  in  (c)  and  with  no  relation  to  corresponding  analyst  phase  notations.  The 
wavefield  counterparts  of  the  P-  and  S-intensity  wavelets  are  given  in  eq.8. 

3.2  Event  location  exercise  using  Italian,  German  and  Norwegian  data 

We  requested  short  period  (z-component)  STA  transformed  data  from  3  local  networks  and  in 
this  respect  as  implied  above  did  not  bother  about  instrument  responses,  local  travel  time  tables, 
and  crustal  models.  In  case  of  the  German  data  we  only  had  at  hand  station  coordinates  and  the 
corresponding  recordings.  The  various  analysis  steps  were  as  follow: 

1.  Band  pass  filtering;  2-4  Hz  pass  band  for  SNR  enhancement.  The  initial  STA-transforms  of 
waveforms  was  (sampling  rate  0.5  sec,  no  overlap):  STA2  =  Ylk  ‘ft 

2.  Derive  model  for  P-  and  S-  mode  intensity  as  a  function  of  distance.  Analysis  step  (1  and 
2)  is  illustrated  in  Figures  4  and  5.  Initially  we  used  Fig.  5a  for  Germany. 

3.  Estimating  minimum  of  entropy  of  source  image  contrast  in  time  (see  Fig.  6  for  visualization) 

4.  Knowing  the  minimum  of  entropy  time  we  searched  the  grid  for  the  most  probable  location 
-  our  choice  of  epicenter  locations  (Fig.  7).  In  Table  1  we  list  reference  epicenter  locations 
(bulletins)  and  relative  positions  of  our  automatic  locations. 

The  following  comments  apply;  our  location  scheme  is  tied  to  epicenter  determinations  from 
which  we  could  accurately  compute  origin  time  for  a  preset  focal  depth  and  local  travel  time  tables. 
Event  and  magnitudes  are  easily  derived  from  the  STA-parameters  using  the  Mendi  and  Husebye 
(1994)  method. 
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Figure  5:  P-  and  S-intensity  wavelet  travel  times  for  Italian,  Norwegian  (a)  and  German  (b) 
network  recordings.  The  corresponding  wavelet  velocities  are  6.3  km/s  and  3.5  km/s  for  (a)  and 
5.9  km/s  and  3.4  km/s  respectively  for  (b). 


Figure  6:  Entropy  of  source  image  contrast  as  a  function  of  time  (eq.15)  for  the  Italian  event 
of  11/02/95  (Table  1).  Note  the  clear  functional  minimum  at  17  sec  being  typical  of  real  event 
detections.  The  corresponding  spatial  distribution  of  the  source  image  function  over  the  spatial 
sampling  area  is  shown  in  Fig.7  . 
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Figure  7:  The  spatial  distribution  of  the  source  image  function  for  the  Italian  event  of  11/02/95 
at  time  17  sec  coinciding  with  the  entropy  minimum  in  Fig.  6.  The  spatial  position  of  the 
functioned  maximum  in  eq.(14)  is  taken  as  the  epicenter  coordinates  (see  listings  in  Table  1,  where 
comparisons  are  made  with  bulletin  locations). 


Date 

Bulletin  data 

Deviation 

Date 

Bulletin  data 

Deviation 

d/m/y 

LatN  LonE 

A  Lai 

A  Lon 

d/m/y/ 

LatN 

LonE 

A  Lat 

A  Lon 

GERMANY 

ITALY 

30/07/94 

51.45 

11.43 

0.06 

0.02 

11/12/94 

37.81 

15.61 

-0.01 

-0.27 

17/09/94 

51.45 

10.50 

0.08 

0.13 

02/01/95 

40.82 

15.34 

-0.15 

mssi 

04/12/94 

50.25 

12.44 

0.23 

Bj&JLJI 

04/01/95 

44.87 

7.36 

-0.11 

-0.26 

16/01/95 

50.83 

6.35 

-0.05 

0.11 

11/02/95 

44.28 

11.50 

0.01 

-0.16 

48.18 

8.75 

-0.11 

-0.13 

18/02/95 

39.14 

16.45 

0.11 

-0.23 

10/03/95 

50.79 

6.35 

0.01 

0.04 

26/02/95 

44.49 

12.08 

-0.05 

0.12 

24/03/95 

47.79 

8.83 

-0.36 

-0.64 

NORWAY 

30/03/95 

51.46 

10.40 

-0.29 

0.39 

26/08/93 

59.05 

5.77 

0.02 

0.24 

51.49 

16.20 

0.18 

-0.16 

26/08/93 

61.11 

4.07 

-0.10 

0.18 

07/06/95 

51.48 

16.20 

0.22 

-0.12 

13/09/93 

66.34 

5.67 

-0.34 

0.68 

Table  1:  Comparisons  between  event  locations  as  reported  in  local  bulletins  and  those  obtained 
using  our  automatic  location  scheme  -  differences  are  mostly  small.  The  few  exceptions  in  this 
regard  are  events  located  outside  or  on  the  periphery  of  the  German  (ca  8  stations  used)  and 
Norwegian  (ca  6  stations  used)  network  areas  where  the  epicenter  resolution  is  relatively  poor  for 
any  location  method,  (ca  20  stations  used  from  Italian  network  area) 


398 


3.3  Research  accomplished 


An  automatic  and  robust  method  for  event  location  using  local  network  data  has  been  developed. 
It  is  easily  adapted  to  any  local  network  area  -  no  need  for  overly  detailed  structural  information. 
It  is  economical  on  two  accounts;  (i)  In  situ  sampling  rate  is  0.5  sec  or  lower  so  minimal  data 
transfer  costs  (ii)  Analyst  interference  with  finalized  event  solutions  is  likely  to  be  minimal  since 
we  bypass  conventional  not  so  robust  analysis  steps  like  signal  detection,  phase  identification  and 
phase  association.  Event  detectability,  not  excessively  tested,  is  likely  to  be  good  as  the  method 
’focuses’  on  maximum  P-and  S-wave-  mode  intensities  and  not  often  weak  phase  onsets  (Fig.  4) 

3.4  Recommendations  and  future  plans 

The  main  recommendation  is  that  the  information  potential  of  local  seismograph  networks  in  the 
context  of  CTBT  monitoring  should  be  exploited  to  an  extent  far  greater  than  currently  under 
consideration.  The  novel  method  of  event  location  presented  provides  one  avenue  for  achieving  a 
better,  faster  and  more  effective  monitoring  environment.  Note,  our  approach  is  based  on  existing 
network  instrumentation  and  perhaps  most  important  does  not  incure  excessive  data  transmission 
costs. 

Options  for  further  method  refinements  are  by  no  means  exhausted;  future  research  topics 
are  better  handling  of  short  distance  records  (P/S  interference/overlap),  interfering  events  and 
naturally  method  adaptation  to  a  real  time  operational  environment.  Future  plans  also  include 
efforts  to  extend  methodology  to  teleseismic  recordings  and  most  challenging  if  source  type  identi¬ 
fication  estimates  can  be  included  for  earthquake  magnitudes  above  3  and  naturally  for  the  many 
stationary  mining  and  quarry  sites. 
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